library(knitr)
## Warning: package 'knitr' was built under R version 4.1.3
opts_chunk$set(tidy.opts = list(width.cutoff = 60))
opts_knit$set(global.par = TRUE)

Reference data

Load Zambia’s ministry of agriculture maize forecasts data as reference.

Load data

## [1] "Province"      "District"      "planted_ha"    "harvested_ha" 
## [5] "production_MT" "yield_MT_ha"   "Year"

Convert District and province names to upper case.

##  [1] "CENTRAL"       "COPPERBELT"    "EASTERN"       "LUAPULA"      
##  [5] "LUSAKA"        "MUCHINGA"      "NORTH-WESTERN" "NORTHERN"     
##  [9] "SOUTHERN"      "WESTERN"
Districts
x
CHADIZA
CHAMA
CHAVUMA
CHIBOMBO
CHIENGE
CHILILABOMBWE

Visualize the crop forecasting data from MoA per year.

boxplot(yield_MT_ha~Year, data=ref, col=rainbow(length(unique(ref$Year))), xlab="Year", ylab = "Yield (MT/ha)", main="Zambia MoA Annual Forecasts.")

Visualize the crop forecasting data from MoA per Province.

Aggregate MoA forecasts per district.

RHEAS simulated yields

Load and aggregate RHEAS simulated Leaf Area Index (LAI), Water stress and Grain Weight Average Dry (GWAD) across different ensembles. Extract year from dates (we will use harvest year).

Unit Production forecast and RHEAS metrics aggregagation

Aggregate RHEAS production forecasts and metrics with respect to Districts maize growing calendar.

The maize growing season in Zambia starts from October to end of June. So we will aggregate the metrics and forecast with this condition using the function RH_metrics.

Convert RHEAS yields from kg/ha to MT/ha.

Validation

Compare RHEAS with MoA forecasts.

Merge the MoA forecasts with RHEAS ones by year and District.

Validation

RMSE

We can use the Root Mean Square Error (RMSE) and mean absolute percentage error (MAPE) to evaluate the models accuracy. RMSE is given as:

\[ \text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^n y-\widehat{y}}, \] where \(\widehat{y}\) and \(y\) are predicted yields and observed yields respectively while n is the number of fitted points.

Compute RMSE.

## RMSE  0.9992922

RRMSE

What about Relative root mean square error (RRMSE):

\[ \text{RRMSE} = \frac{\sqrt{\frac{1}{n} \sum_{i=1}^n y-\widehat{y}}}{\frac{1}{n} \sum_{i=1}^n y}, \]

Compute RRMSE.

## RRMSE  46.20046

According to Li et al. 2013, the performance of the model is excellent when RRMSE < 10%; good if 10% < RRMSE < 20%; fair if 20% < RRMSE < 30%; and poor if RRMSE ≥ 30%.

Total production

Total production P is estimated as a product of unit production u (in MT/ha) and area under maize a per district, that is, $$

p = u a. $$

Merge maize area statistics per district with unit production forecast and compute production in MT/ha.

rh <- merge(rh, a[,c("District","Area_ha")], by="District")
rh$Production <- rh$gwad * rh$Area_ha

Visualization

Add shapefile for visualization.

Check and format District names to be consistent in both the RHEAS and administrative boundaries.

## character(0)
## character(0)

Merge RHEAS and Admin data.

Visualize RHEAS predicted yields spatially.